#ifndef CUDA_FULL_WEIGHTING_C
#define CUDA_FULL_WEIGHTING_C

#include "utility.h"
#include "CUDA_Full_Weighting.h"
#include <math.h>

const int BLOCK_SIZE = 16;
const int R = 1;
__constant__ float weights[3][3];

__global__ void full_weighting_kernel(const float * d_fine, float * d_coarse, int dim_size, int coarse_dim_size)
{
	//Starting at 1 - coarse indices
	int coarse_x = blockIdx.x * blockDim.x + threadIdx.x + 1;
	int coarse_y = blockIdx.y * blockDim.y + threadIdx.y + 1;
	//Starting at 0 - fine indices
	int fine_x = 2 * coarse_x - 1;
	int fine_y = 2 * coarse_y - 1;

	if(coarse_x <= coarse_dim_size && coarse_y <= coarse_dim_size)
	{
		//Compute full weighting
		float value = 0;
		int offset_x, offset_y;
		for(offset_x = -R; offset_x<= R; offset_x++)
		{
			for(offset_y = -R; offset_y<= R; offset_y++)
			{
				value += weights[ offset_x + R ][ offset_y + R ] * d_fine[ (fine_y + offset_y) * dim_size + (fine_x+offset_x) ];
			}
		}
		d_coarse[(coarse_y-1) * coarse_dim_size + (coarse_x-1)] = value;
	}
}

__global__ void full_weighting_kernel_shmem(const float * d_fine, float * d_coarse, int dim_size, int coarse_dim_size)
{
	__shared__ float vec_cache[2*(BLOCK_SIZE + 2*R)][2*(BLOCK_SIZE + 2*R)];
	int coarse_x = blockIdx.x * blockDim.x + threadIdx.x + 1;
	int coarse_y = blockIdx.y * blockDim.y + threadIdx.y + 1;
	//Starting at 0 - fine indices
	int fine_x = 2 * coarse_x - 1;
	int fine_y = 2 * coarse_y - 1;
	if(coarse_x <= coarse_dim_size && coarse_y <= coarse_dim_size)
	{
		int offset_x, offset_y;
		//load values into shared memory

		//full weight computation
		float value = 0;
		for(offset_x = -R; offset_x<= R; offset_x++)
		{
			for(offset_y = -R; offset_y<= R; offset_y++)
			{
				value += weights[ offset_x + R ][ offset_y + R ] * vec_cache[2*threadIdx.x + offset_x + R][2*threadIdx.y + offset_y + R];
			}
		}
		d_coarse[(coarse_y-1) * coarse_dim_size + (coarse_x-1)] = value;
	}
}

void CUDA_Full_Weighting::operator()(const CUDA_vec<float> & a, CUDA_vec<float> & b)
{
	int vec_dim_size = (int)sqrt(a.getSize());
	int coarse_dim = (vec_dim_size+1) / 2 - 1;

	int num_blocks_dim = (int)ceil(((float)coarse_dim) / BLOCK_SIZE);
	dim3 threads_per_block(BLOCK_SIZE,BLOCK_SIZE);
	dim3 blocks_per_grid(num_blocks_dim, num_blocks_dim);
	full_weighting_kernel<<<blocks_per_grid, threads_per_block>>>(a.raw_pointer(),b.raw_pointer(), vec_dim_size, coarse_dim);
	cudaThreadSynchronize();

}



CUDA_Full_Weighting::CUDA_Full_Weighting()
{
	float host_weights[3][3] = { {1/4.0, 1/2.0, 1/4.0}, {1/2.0, 1/1.0, 1/2.0}, {1/4.0, 1/2.0, 1/4.0} };
	cudaMemcpyToSymbol(weights, host_weights, sizeof(host_weights));
}
#endif
